Seattle Pets top of page

Author: Jessica Marx

Date: 01 August 2019

# load packages, install if needed
packages = c(
      "dplyr"
    , "ggplot2"
    , "formattable"
    , "plotly"
    , "RColorBrewer"
    , "scales"
    , "stringr"
    , "tidyr"
    , "ElmeR"
    , "RJDBC"
    , "kableExtra"
    , "wesanderson"
    , "reshape2"
    , "rtweet"
    , "tidytext"
    , "lubridate"
    , "wordcloud"
    , "ggpubr"
    , "ggthemes"
    , "knitrBootstrap"
    , "DT"
    , "MatchIt"
    , "beyonce"
    , "UpSetR"
    , "gganimate"
    , "wordcloud2"
    , "widyr"
    , "ggraph"
    , "igraph"
    , "aod"
    , "corrplot"
    , "ROCR"
    , "InformationValue"
    , "car"
    , "glmnet"
    , "caret"
    , "kernlab"
    , "pdp"
    , "rpart.plot"
    , "rpart"
    , "e1071"
    )
package.check <- lapply(packages, FUN = function(x) {
  if (!require(x, character.only = TRUE)) {
    install.packages(x, dependencies = TRUE)
    library(x, character.only = TRUE)
  }
})
options(scipen= 999)
theme_set(theme_minimal(base_size = 9, base_family = "Roboto"))

#functions!
#round 
round.to <- function(x, b) {
  round(x/b)*b
}
#odds to probability
odds.to.prob <- function(odds) {
  odds/(1 + odds) 
}
#log odds to probability 
logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}
#convert to a range
range01 <- function(x){
  (x-min(x))/(max(x)-min(x))
}
#function to get vector of color values from RColorBrewer
get_hex_values <- function(pal) {
    brewer.pal(brewer.pal.info[pal, "maxcolors"], pal)
}
paired_cols <- get_hex_values(pal = "Paired")

Top Dog Breeds

a = seattle_pet_licenses %>% 
  filter(species == "Dog") %>% 
  select(primary_breed, license_issue_date) %>% 
  mutate(license_issue_date = as.Date(license_issue_date),
         year = lubridate::year(license_issue_date)) %>% 
  group_by(primary_breed, year) %>% 
  count() %>% 
  arrange(year, desc(n)) %>% 
  ungroup() %>% 
  group_by(year) %>% 
  mutate(rank = row_number()) %>% 
  filter(rank <= 5, year >= 2010) %>%
  mutate(rev_rank = rank(-rank)) %>% 
  ungroup() 
colourCount = length(unique(a$primary_breed))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
ggplot(a, aes(x = year, y = reorder(rev_rank, rev_rank), fill = primary_breed)) +
  geom_bar(stat = "identity", position = "stack", color = "white") + 
  geom_text(stat = "identity", position = "stack", aes(label = rank), color = "black", size = 3, vjust = 1, hjust = 3) + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_x_continuous(breaks = seq(2010, 2016, 1)) + 
  coord_flip() +
  ylab("rank") + 
  scale_fill_manual(values = getPalette(colourCount), name = "Breed") +  
  ggtitle(label = "Top Five Seattle Dog Breeds by Year")

Breed Acceleration

a = seattle_pet_licenses %>% 
  filter(species == "Dog") %>% 
  select(primary_breed, license_issue_date) %>% 
  mutate(license_issue_date = as.Date(license_issue_date),
         year = lubridate::year(license_issue_date)) %>% 
  filter(year >= 2015) %>% 
  group_by(primary_breed, year) %>% 
  count() %>%
  spread(year, n) %>% 
  drop_na() %>% 
  mutate(pct_increase = `2016`/`2015`-1) %>% 
  arrange(desc(pct_increase)) %>% 
  head(20) %>% 
  ungroup()
b = seattle_pet_licenses %>% 
  filter(species == "Dog") %>% 
  select(secondary_breed, license_issue_date) %>% 
  mutate(license_issue_date = as.Date(license_issue_date),
         year = lubridate::year(license_issue_date)) %>% 
  filter(year >= 2015) %>% 
  group_by(secondary_breed, year) %>% 
  count() %>%
  spread(year, n) %>% 
  drop_na() %>% 
  mutate(pct_increase = `2016`/`2015`-1) %>% 
  arrange(desc(pct_increase)) %>% 
  head(20) %>% 
  ungroup()
breeds = a %>% 
  select(primary_breed, pct_increase) %>% 
  rename("breed" = primary_breed) %>% 
  mutate(breed_level = "Primary Breed") %>% 
  rbind(
    b %>% 
  select(secondary_breed, pct_increase) %>% 
  rename("breed" = secondary_breed) %>% 
  mutate(breed_level = "Secondary Breed")
  ) %>% 
  group_by(breed_level) %>% 
  arrange(breed_level, desc(pct_increase)) %>% 
  mutate(rank = row_number()) %>% 
  ungroup()
  
colourCount = length(unique(breeds$breed))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
breeds %>% 
  ggplot(aes(x = reorder(breed, -rank), y = pct_increase, fill = breed)) +
  geom_bar(stat = "identity") + 
  scale_fill_manual(values = getPalette(colourCount), name = "Breed") +  
  scale_y_continuous(labels = percent_format()) +
  ylab("Percent Increase") + 
  xlab("Breed") + 
  ggtitle(label = "Year Over Year (YOY) Increase in Breeds") + 
    coord_flip() + 
    theme(legend.position = "none") + 
  facet_wrap(~breed_level, scales = "free_y")

Dog Names

a = seattle_pet_licenses %>% 
  filter(species == "Dog",
         !is.na(animal_s_name)) %>% 
  select(animal_s_name, license_issue_date) %>% 
  mutate(animal_s_name = if_else(is.na(animal_s_name), "No Name", animal_s_name) ,
         license_issue_date = as.Date(license_issue_date),
         year = lubridate::year(license_issue_date)) %>% 
  group_by(animal_s_name, year) %>% 
  count() %>% 
  arrange(year, desc(n)) %>% 
  ungroup() %>% 
  group_by(year) %>% 
  mutate(rank = row_number()) %>% 
  filter(rank <= 5, year >= 2010) %>%
  mutate(rev_rank = rank(-rank)) %>% 
  ungroup() 
colourCount = length(unique(a$animal_s_name))
getPalette = colorRampPalette(brewer.pal(12, "Paired"))
ggplot(a, aes(x = year, y = reorder(rev_rank, rev_rank), fill = animal_s_name)) +
  geom_bar(stat = "identity", position = "stack", color = "white") + 
  scale_fill_manual(values = getPalette(colourCount), name = "Name") +  
  geom_text(stat = "identity", position = "stack", aes(label = rank), color = "black", size = 3, vjust = 1, hjust = 3) + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  scale_x_continuous(breaks = seq(2010, 2016, 1)) + 
  coord_flip() +
  ylab("rank") + 
  
  ggtitle(label = "Top Five Seattle Dog Names by Year")

Breed Combos

Breed Heatmap

Which primary and secondary breek combinations are most popular?

top_breeds = rbind(
(seattle_pet_licenses %>% 
  filter(species == "Dog") %>% 
  select(primary_breed) %>% 
  group_by_all() %>% 
  count() %>% 
  arrange(desc(n)) %>% 
  head(20) %>% 
  rename("breed" = primary_breed)),
(seattle_pet_licenses %>% 
  filter(species == "Dog") %>% 
  select(secondary_breed) %>% 
  group_by_all() %>% 
  count() %>% 
  arrange(desc(n)) %>% 
  head(20) %>% 
  rename("breed" = secondary_breed)
)) %>% 
  group_by(breed) %>% 
  summarise(n = sum(n)) %>% 
  arrange(desc(n)) %>% 
  ungroup() %>% 
  mutate(breed = if_else(is.na(breed), "Unknown", breed))
plot = seattle_pet_licenses %>% 
  filter(species == "Dog", 
         primary_breed %in% top_breeds$breed
         ,secondary_breed %in% top_breeds$breed
         ) %>% 
  select(primary_breed, secondary_breed) %>% 
  mutate(secondary_breed = if_else(is.na(secondary_breed), primary_breed, secondary_breed)) %>% 
  group_by_all() %>% 
  count() %>% 
  arrange(desc(n)) %>% 
  #head(50) %>% 
  mutate(log = log(n + 1)) %>% 
  ggplot(
    aes(x = primary_breed,
        y = secondary_breed, 
        fill = log,
        text = paste(
          "Primary Breed:", primary_breed,
          "<br>Secondary Breed:", secondary_breed,
          "<br>Dogs:", comma(n)
        )
      )
    ) + 
  geom_tile() +
  theme(axis.text.x = element_text(angle = -30, hjust = 0)) + 
  scale_fill_gradientn(colours = blues9,
    na.value = "white",
    breaks=c(0,7),
    # # # #breaks = c(400, 700000),
    labels=c("0","1K"),
    limits=c(0,7),
    # # # #limits = c(400, 700000),
    name = "Dogs") +
  ylab("Secondary Breed") +
  xlab("Primary Breed") 
ggplotly(plot, tooltip = "text") %>% 
  layout(autosize = F, width = 800, height = 600)

Upset Plot (like a Venn Diagram, but better)

---

#This is the true code for this plot. Unfortunately, this package has a bug -- one that prints a giant blank space before the plot, which makes it less than optimal for publishing. My hack around this is to reproduce the plot using this code and save the image in my R Project file. From there I can just pull it in with html. A hack, yes, but one that works until the bug is fixed (according to their GitHub page, it's a WIP!). 
# upset_dog = rbind(
# (seattle_pet_licenses %>% 
#   filter(species == "Dog") %>% 
#    mutate(primary_breed = if_else(is.na(primary_breed), "Unknown", primary_breed)) %>% 
#   select(primary_breed, license_number) %>% 
#   rename("breed" = primary_breed)),
# (seattle_pet_licenses %>% 
#   filter(species == "Dog") %>% 
#    mutate(secondary_breed = if_else(is.na(secondary_breed), "Unknown", secondary_breed)) %>%  
#   select(secondary_breed, license_number) %>% 
#   rename("breed" = secondary_breed)
# )) %>% 
#   group_by(license_number, breed) %>%
#   count() %>% 
#   mutate(value = 1) %>%
#   select(-license_number, -n) %>% 
#   spread(breed, value, fill = 0) %>%
#   as.data.frame()
# 
# upset(upset_dog, 
#       order.by = "freq",
#       nsets = 15, 
#       mainbar.y.label = "Dog Breed Intersections", 
#       sets.x.label = "Total Number of Breed in Sample",
#       matrix.color = "#3288bd",
#       main.bar.color = "#f46d43",
#       sets.bar.color = "#3288bd", 
#       point.size = 2.5, line.size = 1, 
#       nintersects = 15,
#       shade.color = "#9e0142",
#       text.scale = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
#       group.by = "degree")

Breed Network

Breed Characteristics

What are the words used to describe different dog breeds? All breed descriptions were obtained through the Pet Breeds Characteristics public dataset.

We will score the top 20 breeds using the AFINN lexicon, which scores words by sentiment.

Before doing so, let’s inspect the data and make sure that we agree with the way the Lexicon is classifying words given our specific context.

AFINN <- get_sentiments("afinn")
pal <- rev(beyonce_palette(22, 100, type = "continuous"))
plot_c = barsentiment %>%
    dplyr::count(sentiment, word, n=n) %>%
    dplyr::ungroup() 
dog_breed_characteristics %>% 
  select(BreedName, Group1, Group2, Temperment) %>% 
  unnest_tokens(word, Temperment) %>% 
  inner_join(AFINN) %>% 
  select(word, score) %>% 
  unique() %>% 
  arrange(desc(score)) %>% 
  mutate(Sentiment = if_else(score > 0, "Positive", "Negative")) %>% 
  ggplot(aes(x = reorder(word, score), y = score, fill = Sentiment)) +
  geom_bar(stat = "identity") + 
  scale_fill_brewer(palette = "Set1") + 
  coord_flip() + 
  xlab("Word") + 
  ylab("Score")

A dog being alert does not seem like it would be a bad thing - let’s change that from a -1 to a +1.

Some dog breeds are described by more words than others. A histogram helps us visualize and normalize the data.

Now we’re ready to score our dogs by breed.

Top 20

dog_scores %>% 
  arrange(desc(total_score)) %>% 
  head(20) %>% 
  ggplot(aes(x = reorder(BreedName, total_score), y = total_score, fill = total_score)) +
  geom_bar(stat = "identity") + 
  scale_fill_gradientn(
    colours = pal, name = "Score") +
  coord_flip() + 
  xlab("Breed") + 
  ylab("Total Word Score")

Are dogs with more words scored higher? We can normalize by taking the mean word score.

dog_scores %>% 
  arrange(desc(mean_score), BreedName) %>%  
  head(20) %>% 
  ggplot(aes(x = reorder(BreedName, mean_score), y = mean_score, fill = total_score)) +
  geom_bar(stat = "identity") + 
  scale_fill_gradientn(
    colours = pal, name = "Total Score") +
  coord_flip() + 
  xlab("Breed") + 
  ylab("Mean Word Score")

Top 20 Breed Combos

By Group

Let’s hear it for the the Hounds!

group_scores = dog_words %>% 
  filter(!Group1 == "Southern") %>% 
  select(Group1, word) %>% 
  unique() %>% 
  inner_join(AFINN) %>% 
  group_by(Group1) %>% 
  summarise(avg_score = mean(score),
            score = sum(score)) %>% 
  drop_na()
group_scores %>% 
  arrange(desc(avg_score)) %>% 
  #head(20) %>% 
  ggplot(aes(x = reorder(Group1, avg_score), y = avg_score, fill = score)) +
  geom_bar(stat = "identity") + 
  scale_fill_gradientn(
    colours = pal, name = "Total Score") +
  coord_flip() + 
  xlab("Breed Group") + 
  ylab("Mean Score")

Breed Group 1 + Group 2

group_both = dog_words %>% 
  filter(!Group1 == "Southern") %>% 
  inner_join(AFINN) %>% 
  mutate(group_combo = paste(Group1, Group2, sep = " + ")) %>% 
  group_by(group_combo) %>% 
  summarise(
    avg_score = mean(score), 
    score = sum(score)
    ) %>% 
  #summarise(count = n_distinct(word)) %>% 
  drop_na()
group_both %>% 
  arrange(desc(avg_score)) %>% 
  #head(20) %>% 
  ggplot(aes(x = reorder(group_combo, avg_score), y = avg_score, fill = score)) +
  geom_bar(stat = "identity") + 
  scale_fill_gradientn(
    colours = pal, name = "Total Score") +
  coord_flip() + 
  xlab("Breed Group 1 + Breed Group 2") + 
  ylab("Mean Score")

Intelligence Ratings

Does intelligence correspond with price? No.

Pets by Zip YOY

seattle_pet_licenses %>% 
    drop_na() %>% 
  select(license_issue_date
         , species
         , zip_code
         , license_number) %>% 
  mutate(license_issue_date = as.Date(license_issue_date),
         year = lubridate::year(license_issue_date),
         zip_code = (as.factor(zip_code))
         ) %>%
  group_by(year, zip_code, species) %>% 
  summarise(new_pets = n_distinct(license_number)) %>% 
  ungroup() %>% 
  plot_ly(
  x = ~zip_code
  , y = ~new_pets
  , color = ~species
  , type = "bar"
  , legend_group = ~species
  , frame = ~year
  ) %>% 
  layout(
    barmode = "stack"
    , xaxis = list(title = "Zip Code")
    , yaxis = list(title = "New Pets")
  )